rm(list = ls())
library(cowplot)
library(rootSolve)
library(latex2exp)
alpha <- 0.36 # 中间产品参数
theta <- 1.5 # 消费参数
rho <- 0.01 # 贴现率

A1 <- 1
L1 <- 1
L2 <- 1
A2 <- 0.5

# 影响渠道的正反通过sig1正负来设置
sig1 <- 0.4 # c1的研发成本参数，对应epu1
sig2 <- 0.6 # c2的研发成本参数
sig3 <- 0.4 # c2的研发成本参数,对应epu2
sig4 <- 0.7
phi1 <- 1 # 研发成本参数
phi2 <- 1

Fn <- function(x, phi1, phi2, epu1, sig1, sig2, sig3, L1, L2, A1, A2,alpha,theta,rho){
  # 中间变量计算
  epu2 <- 0.5 * epu1
  pi1 <- L1*A1^(1/(1-alpha))*((1-alpha)/alpha)*alpha^(2/(1-alpha)) # 利润
  pi2 <- L2*A2^(1/(1-alpha))*((1-alpha)/alpha)*alpha^(2/(1-alpha)) # 利润
  eta1 <- phi1 * epu1^sig1 * x[1]^sig2 # 研发成本
  eta2 <- phi2 * epu2^sig3 * (x[3]/x[1])^sig4 # 研发成本
  # 方程组
  # N
  f1 <- 1/eta1 * (A1^(1/(1-alpha)) * L1 * alpha^(2*alpha/(1-alpha)) * (1-alpha^2) - x[2]/x[1])
  # C
  f2 <- 1/theta * (pi1/eta1-rho)
  f3 <- 1/eta2 * (A2^(1/(1-alpha)) * L2 * alpha^(2*alpha/(1-alpha)) * (1-alpha^2) - x[4]/x[3])
  f4 <- 1/theta * (pi2/eta2-rho)
  return(c(f1 = f1, f2 = f2, f3 = f3, f4 = f4))
}

# 设置初值
x <- c(N1 = 1, cons1 = 0.5, N2 = 0.9, cons2 = 0.3)

DrawSimu <- function(sig1 = 0.4,sig3 = 0.4, tt = ''){
  picdata <- data.frame(epu1 = seq(1,4,0.1), N2 = NA)
  for (i in 1:nrow(picdata)) {
    rlt <- multiroot(f = Fn, start = x, phi1=phi1, phi2=phi2, epu1= picdata$epu1[i], sig1=sig1,
                     sig2= sig2, sig3=sig3, L1=L1, L2=L2, A1=A1, A2=A2,alpha=alpha,theta=theta,rho=rho)
    picdata$N2[i] <- rlt$root[3]
  }
  library(ggplot2)
  p <- ggplot(picdata, aes(x = epu1, y = N2)) + geom_line() +
    labs(y = TeX('$N_2$'), x= TeX('$EPU_1$'), title = tt) + theme_bw()+
    theme(plot.title = element_text(hjust = 0.5))
  return(p)
}

p1 <- DrawSimu(sig1 = 0.4,sig3 = 0.3, tt = TeX('$(a)\\; \\sigma_2 = 0.4,\\; \\sigma_4 = 0.3$'))
p2 <- DrawSimu(sig1 = -0.4,sig3 = 0.3,tt = TeX('$(b)\\;\\sigma_2 = -0.4,\\;\\sigma_4 = 0.3$'))
p22 <- DrawSimu(sig1 = -0.4,sig3 = 0.6, tt = TeX('$(c)\\;\\sigma_4 = -0.4,\\;\\sigma_4 = 0.6$'))
p5 <- plot_grid(p1,p2,p22,ncol = 1)
# ggsave('data-raw/simu1.wmf')

p3 <- DrawSimu(sig3 = 0.5, sig1 = 0.4,tt = TeX('$(d)\\;\\sigma_2 = 0.4, \\;\\sigma_4 = 0.5$'))
p4 <- DrawSimu(sig3 = -0.5, sig1 = 0.4,tt = TeX('$(e)\\;\\sigma_2 = 0.4, \\;\\sigma_4 = -0.5$'))
p42 <- DrawSimu(sig1 = -0.5,sig3 = 0.6,tt = TeX('$(f)\\;\\sigma_2 = 0.6, \\;\\sigma_4 = -0.5$'))
p6 <- plot_grid(p3,p4,p42,ncol = 1)
plot_grid(p5,p6,ncol = 2)
# ggsave('./data-raw/simuall.wmf')
